Packages

pacman::p_load(
  dplyr,
  tidyverse,
  readxl,
  vegan)

# setwd("D:PhD/01_Data/03_Vitro/02_Identification of glycan utilizing microbe") # Windows
setwd("/Volumes/Yiming_Wang/PhD/01_Data/03_Vitro/02_Identification of glycan utilizing microbe") # Mac
# setwd("/Users/godric/Desktop/Fut2 review comments/Revision") #Macbook pro

df_all <- read_excel("Glycan utilization_L6_taxa formatted.xlsx") %>% 
  mutate(Genotype=factor(Genotype,levels = c("WT","KO"))) %>% 
  mutate(Sex = factor(Sex, levels = c("Male", "Female"))) %>% 
  mutate(Genotype_Sex = factor(Genotype_Sex,levels = c("WT Male","WT Female","KO Male","KO Female"))) %>% 
  filter(Genotype == "WT") %>% 
  mutate(Sample = row_number()) %>% 
  mutate(Treatment = ifelse(Treatment == "No glycan", "mBasal", "mBasal + 2'-FL")) %>% 
  mutate(Treatment_Sex= paste0(Treatment,"_", Sex)) %>% 
  relocate(Sample, Treatment_Sex)

# Change 1 to 01
df_all$Sample <- sprintf("%02d", as.numeric(df_all$Sample)) 
df_all<- df_all %>% 
  mutate(Sample = factor(paste0("Sample", " ",df_all$Sample)))

# Subgroup by variables
#30436d WT Male
#9ea5c2 WT Female
#654e3a KO Male
#bcaba6 KO Female

#3C5488B2 WT
#7E6148B2 KO
#D98594FF Female
#94C5CCFF Male

Alpha diversity_visualization and analysis

Plot

# load package
pacman::p_load(
  ggplot2,ggsci, ggthemes,
)

set.seed(123)

# Level orders
level_order_Treatment <- c("mBasal + 2'-FL", "mBasal")

# My color panel
mypal = pal_npg("nrc", alpha = 0.7)(9)
mypal
## [1] "#E64B35B2" "#4DBBD5B2" "#00A087B2" "#3C5488B2" "#F39B7FB2" "#8491B4B2"
## [7] "#91D1C2B2" "#DC0000B2" "#7E6148B2"
library("scales")
show_col(mypal) #  #3C5488B2(WT) and  #E64B35B2 (KO)

# Treatment
## FDP Diversity
Treatment_FDP <- ggplot (data=df_all, mapping = aes(y=FPD_diversity, x=factor(Treatment, level=level_order_Treatment))) +
  geom_boxplot(outlier.shape = NA)+
  theme_bw() +
  geom_jitter(aes(colour=Treatment),size=2,alpha=0.9,
              position=position_jitterdodge(jitter.width = 0.5,
                                            jitter.height = 0, 
                                            dodge.width = 0.8))+ 
  labs(x="", y="Faith's phylogenetic diversity", caption ="", # x-,y-axis name
       title = "")+
  theme(legend.position = "none", # right,left, bottom
        legend.title = element_blank(), # can remove this
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 5)),
        axis.title.y=element_text(margin = margin(r = 5), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_y_continuous(limits = c(2,10),breaks = seq(2,10,2))+
  scale_colour_manual(values=c("#b86877","#70a1a7")) + 
  annotate("text",x=1.5,y=10,
           label="p = 0.0094",hjust=0.5, size=5.5)


## Shannon Diversity
Treatment_shannon <- ggplot (data=df_all, mapping = aes(y=Shannon_diversity, x=factor(Treatment, level=level_order_Treatment))) +
  geom_boxplot(outlier.shape = NA)+
  theme_bw() +
  geom_jitter(aes(colour=Treatment),size=2,alpha=0.9,
              position=position_jitterdodge(jitter.width = 0.5,
                                            jitter.height = 0, 
                                            dodge.width = 0.8))+ 
  labs(x="", y="Shannon diversity", caption ="", # x-,y-axis name
       title = "")+
  theme(legend.position = "none", # right,left, bottom
        legend.title = element_blank(), # can remove this
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 5)),
        axis.title.y=element_text(margin = margin(r = 5), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_y_continuous(limits = c(3.5,5),breaks = seq(3.5,5,0.3))+
  scale_colour_manual(values=c("#b86877","#70a1a7")) + 
  annotate("text",x=1.5,y=5,
           label="p = 0.0020",hjust=0.5, size=5.5)


## Richness (Chao)
Treatment_chao <- ggplot (data=df_all, mapping = aes(y=Chao_richness, x=factor(Treatment, level=level_order_Treatment))) +
  geom_boxplot(outlier.shape = NA)+
  theme_bw() +
  geom_jitter(aes(colour=Treatment),size=2,alpha=0.9,
              position=position_jitterdodge(jitter.width = 0.5,
                                            jitter.height = 0, 
                                            dodge.width = 0.8))+ 
  labs(x="", y="Chao richness", caption ="", # x-,y-axis name
       title = "")+
  theme(legend.position = "none", # right,left, bottom
        legend.title = element_blank(), # can remove this
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 5)),
        axis.title.y=element_text(margin = margin(r = 5), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_y_continuous(limits = c(30,50),breaks = seq(30,50,5))+
  scale_colour_manual(values=c("#b86877","#70a1a7")) + 
  annotate("text",x=1.5,y=50,
           label="ns",hjust=0.5, size=5.5)



## Pielou evenness
Treatment_pielou <- ggplot (data=df_all, mapping = aes(y=Pielou_evenness, x=factor(Treatment, level=level_order_Treatment))) +
  geom_boxplot(outlier.shape = NA)+
  theme_bw() +
  geom_jitter(aes(colour=Treatment),size=2,alpha=0.9,
              position=position_jitterdodge(jitter.width = 0.5,
                                            jitter.height = 0, 
                                            dodge.width = 0.8))+ 
  labs(x="", y="Pielou evenness", caption ="", # x-,y-axis name
       title = "")+
  theme(legend.position = "none", # right,left, bottom
        legend.title = element_blank(), # can remove this
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 5)),
        axis.title.y=element_text(margin = margin(r = 5), size=16),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_y_continuous(limits = c(0.7,0.9),breaks = seq(0.7,0.9,0.05))+
  scale_colour_manual(values=c("#b86877","#70a1a7")) + 
  annotate("text",x=1.5,y=0.9,
           label="p < 0.0001",hjust=0.5, size=5.5)






# Print
Treatment_FDP

Treatment_shannon

Treatment_chao

Treatment_pielou

Beta diversity_results and visualisation

Datashape

# places that are not satisfied
  ## 1) Permanova results table
  ## 2) Figure A and B did not align because of legend (Change legend size, position to bottom, change P value position)


#library packages
pacman::p_load(
  data.table,
  vegan,
  knitr,
  EcolUtils
)

# Shape data
# Level orders
df_all$Treatment <- factor (df_all$Treatment, levels = c("mBasal + 2'-FL", "mBasal"))
df_all$Treatment_Sex <- factor (df_all$Treatment_Sex, levels = c("mBasal + 2'-FL_Male","mBasal + 2'-FL_Female", "mBasal_Male", "mBasal_Female"))


# subset dateframe
df_abundance <- subset(df_all, select = -c(1:15))
df_metadata <- subset(df_all, select = c(1:15))

Bray-curtis

 # Check the stress and choose the optimal number of dimensions
dist_bray <- vegdist(df_abundance,  method = "bray")


# NMDS ordination_bray curtis
NMDS <- metaMDS(df_abundance, distance = "bray", k = 2)
## Run 0 stress 0.01751042 
## Run 1 stress 0.02006365 
## Run 2 stress 0.01758811 
## ... Procrustes: rmse 0.02845926  max resid 0.08098576 
## Run 3 stress 0.01751045 
## ... Procrustes: rmse 0.0004187265  max resid 0.001625643 
## ... Similar to previous best
## Run 4 stress 0.0190382 
## Run 5 stress 0.01962485 
## Run 6 stress 0.01655679 
## ... New best solution
## ... Procrustes: rmse 0.01742603  max resid 0.04483727 
## Run 7 stress 0.01792821 
## Run 8 stress 0.01704967 
## ... Procrustes: rmse 0.01017602  max resid 0.03491165 
## Run 9 stress 0.01655727 
## ... Procrustes: rmse 0.0002291111  max resid 0.0009842971 
## ... Similar to previous best
## Run 10 stress 0.01790897 
## Run 11 stress 0.01790904 
## Run 12 stress 0.01951685 
## Run 13 stress 0.02049977 
## Run 14 stress 0.01751057 
## Run 15 stress 0.01775668 
## Run 16 stress 0.01879263 
## Run 17 stress 0.02018922 
## Run 18 stress 0.01981532 
## Run 19 stress 0.01792839 
## Run 20 stress 0.01694019 
## ... Procrustes: rmse 0.02756249  max resid 0.08516407 
## *** Solution reached
NMDS$stress
## [1] 0.01655679
# PERMANOVA and Pair-wise PERMANOVA

set.seed(123)
## Genotype and sex are crossed variable, cage are nested in genotype and sex
adonis2(df_abundance ~ Treatment, data=df_metadata, permutations=9999, method="bray")
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## adonis2(formula = df_abundance ~ Treatment, data = df_metadata, permutations = 9999, method = "bray")
##           Df SumOfSqs     R2      F Pr(>F)    
## Treatment  1  0.85111 0.8794 131.25  1e-04 ***
## Residual  18  0.11672 0.1206                  
## Total     19  0.96783 1.0000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(df_abundance ~ Treatment_Sex, data=df_metadata, permutations=9999, method="bray") 
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
## 
## adonis2(formula = df_abundance ~ Treatment_Sex, data = df_metadata, permutations = 9999, method = "bray")
##               Df SumOfSqs      R2      F Pr(>F)    
## Treatment_Sex  3  0.86871 0.89758 46.742  1e-04 ***
## Residual      16  0.09912 0.10242                  
## Total         19  0.96783 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df_permanova_pairwise <- adonis.pair(vegdist(df_all[,16:57],  method = "bray"), df_all$Treatment_Sex, nper = 9999, corr.method = "fdr")
kable(df_permanova_pairwise)
combination SumsOfSqs MeanSqs F.Model R2 P.value P.value.corrected
mBasal + 2’-FL_Male <-> mBasal + 2’-FL_Female 0.0150608 0.0150608 1.9314491 0.1944781 0.1178 0.14136
mBasal + 2’-FL_Male <-> mBasal_Male 0.3907454 0.3907454 52.3473072 0.8674340 0.0103 0.01545
mBasal + 2’-FL_Male <-> mBasal_Female 0.3846677 0.3846677 45.2388265 0.8497337 0.0085 0.01545
mBasal + 2’-FL_Female <-> mBasal_Male 0.4782065 0.4782065 123.0189156 0.9389401 0.0082 0.01545
mBasal + 2’-FL_Female <-> mBasal_Female 0.4662045 0.4662045 94.6449896 0.9220615 0.0071 0.01545
mBasal_Male <-> mBasal_Female 0.0025384 0.0025384 0.5527043 0.0646233 0.7113 0.71130
# Assessment of multivariate homogeneity of groups dispersions
dispersion_bray <- betadisper(dist_bray, group = df_metadata$Treatment)
permutest(dispersion_bray)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df    Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.0014151 0.0014151 1.0824    999  0.356
## Residuals 18 0.0235310 0.0013073
# Figures
#library packages
pacman::p_load(
  ggthemes,
  ggsci
)
## Set seed
set.seed(123)
# options( scipen = 999 )

# Genotype
datascores <- as.data.frame(scores(NMDS))#$sites
scores <- cbind(as.data.frame(datascores), Treatment = df_metadata$Treatment)
centroids <- aggregate(cbind(NMDS1, NMDS2) ~ Treatment, data = scores, FUN = mean)
seg <- merge(scores, setNames(centroids, c('Treatment','oNMDS1','oNMDS2')),
             by = c('Treatment'), sort = FALSE)

composition_Treatment <- ggplot(scores, aes(x = NMDS1, y = NMDS2, color = Treatment)) + 
  geom_segment(data = seg,
             mapping = aes(xend = oNMDS1, yend = oNMDS2),
             size=0.25) + 
  geom_point(data = centroids, size = 4) +                    
  geom_point()+
  coord_fixed()+ 
  theme_bw() +
  labs(title="")+
  theme(legend.position = "bottom",
        legend.title= element_text(size=10, face="plain"),
        legend.text = element_text(size=14),
        legend.key = element_rect(fill = "transparent", colour = "black"),
        axis.title = element_text(face="plain", size = 10),
        axis.title.x=element_text(margin = margin(t = 9), size=14),
        axis.title.y=element_text(margin = margin(r = 5), size=14),
        axis.text.x=element_text(colour="black", face="plain", size=14), 
        axis.text.y=element_text(colour="black", face="plain", size=14),
        plot.title = element_text(size=10, hjust=0.5, face="plain"),
        panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
                                    fill = NA,
                                    size = 1),
        legend.background = element_rect(color="transparent"),
        panel.grid = element_blank(),
        aspect.ratio = 1,
        plot.background = element_rect(fill="transparent", colour = "transparent"))+
  scale_x_continuous(limits = c(-0.5,0.5),breaks = seq(-0.5,0.5,0.25))+
  scale_y_continuous(limits = c(-0.25,0.25),breaks = seq(-0.25,0.25,0.1))+
  scale_color_manual(values=c("#70a1a7", "#b86877"))+
  scale_fill_manual(values=c("#70a1a7", "#b86877"))+
  annotate("text",x=-0.5,y=-0.25, #0.4 if legend position is right
           label="Stress = 0.017",hjust=0, size=4.5)+ #2.2 if legend position is right
  annotate("text",x=-0,y=0.25,
           label=expression(italic("P") < "0.0001"),hjust=0.5, size=4.5)

composition_Treatment

ANOISM: Variance between groups and within groups

# Bray-curtis
dist_bray <- vegdist(df_abundance,  method = "bray")
## Genotype
anosim.Treatment_bray <- with(df_metadata, anosim(dist_bray, Treatment))
summary(anosim.Treatment_bray)
## 
## Call:
## anosim(x = dist_bray, grouping = Treatment) 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.998 
##       Significance: 0.001 
## 
## Permutation: free
## Number of permutations: 999
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0971 0.1675 0.2408 0.3287 
## 
## Dissimilarity ranks between and within classes:
##                0%    25%   50%    75% 100%   N
## Between        87 115.75 140.5 165.25  190 100
## mBasal + 2'-FL  1  29.00  56.0  79.00   93  45
## mBasal          3  21.00  40.0  57.00   77  45
plot(anosim.Treatment_bray)